rm(list = ls())
require(stats4)
require(nlWaldTest)
require("multiwayvcov")
require(msm)
require(stargazer)
require(mfx)

data21=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december6noon_type7clean.csv", header=TRUE)
data22=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december15_type7clean.csv", header=TRUE)
data23=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_january26_type7(clean).csv", header=TRUE)

dataframe21=data.frame(data21)
dataframe22=data.frame(data22)
dataframe23=data.frame(data23)

#some definitions
expframe=rbind(dataframe21,dataframe22,dataframe23)
expframe=subset(expframe,uid>800)
expframe$correct <-(expframe$state==4&(expframe$response==2|expframe$response==7|expframe$response==8|expframe$response==9))|(expframe$state==3&(expframe$response==1|expframe$response==6|expframe$response==3|expframe$response==5))
expframe$incentive <- (expframe$qid==5)*5+(expframe$qid==1)*40+(expframe$qid==6)*70+(expframe$qid==7)*95
accuracy=c(sum((expframe$incentive==5)*expframe$correct)/sum(expframe$incentive==5),sum((expframe$incentive==40)*expframe$correct)/sum(expframe$incentive==40),sum((expframe$incentive==70)*expframe$correct)/sum(expframe$incentive==70),sum((expframe$incentive==95)*expframe$correct)/sum(expframe$incentive==95))
incentivec=c(5,40,70,95)
uidlist=unique(expframe$uid)

#check N
#**
length(uidlist)
#**

#for ease of factor use
expframe$incentive5=(expframe$incentive==5)
expframe$incentive40=(expframe$incentive==40)
expframe$incentive70=(expframe$incentive==70)
expframe$incentive95=(expframe$incentive==95)

#Check NIAS 
expframe$Bresponse <-(expframe$response==2|expframe$response==7|expframe$response==8|expframe$response==9)
expframe$Bstate <- expframe$state==4
expframe$Aresponse <-(expframe$response==1|expframe$response==6|expframe$response==3|expframe$response==5)
expframe$Astate <- expframe$state==3

#NIAS in aggregate

#Base Regression
expframe$incentivedec=expframe$incentive/100
modelNIAS=lm(Aresponse~factor(Astate+incentivedec)+0,data=expframe)
modelNIAC=modelNIAS

#Prob correct vs incentive
fitun=vector('numeric')

#Adjust COV
vcovmodNIAS=vcovHC(modelNIAS, type="HC3")
vcovmodNIAScluster=cluster.vcov(modelNIAS, expframe$uid)


fitun=vector("numeric")
fitun[1]=(modelNIAC$coefficients[5]+1-modelNIAC$coefficients[1])/2
fitun[2]=(modelNIAC$coefficients[6]+1-modelNIAC$coefficients[2])/2
fitun[3]=(modelNIAC$coefficients[7]+1-modelNIAC$coefficients[3])/2
fitun[4]=(modelNIAC$coefficients[8]+1-modelNIAC$coefficients[4])/2


Textstr=vector("character")
Pvalvec=vector("numeric")

#Tests
Textstr[1]="(b[5]-b[1])"
wald1=nlWaldtest(modelNIAS,texts=Textstr,Vcov=vcovmodNIAScluster)
Pvalvec[1]=wald1$p.value

Textstr[1]="(b[6]-b[2])"
wald2=nlWaldtest(modeNIAS,texts=Textstr,Vcov=vcovmodNIAScluster)
Pvalvec[2]=wald1$p.value

Textstr[1]="(b[7]-b[3])"
wald3=nlWaldtest(modelNIAS,texts=Textstr,Vcov=vcovmodNIAScluster)
Pvalvec[3]=wald1$p.value

Textstr[1]="(b[8]-b[4])"
wald4=nlWaldtest(modelNIAS,texts=Textstr,Vcov=vcovmodNIAScluster)
Pvalvec[4]=wald1$p.value

Textstr=vector("character")
#Joint test
Textstr[1]="(b[5]-b[1])"
Textstr[2]="(b[6]-b[2])"
Textstr[3]="(b[7]-b[3])"
Textstr[4]="(b[8]-b[4])"
nlWaldtest(modelNIAS,texts=Textstr,Vcov=vcovmodNIAScluster)

pAAvec=c(modelNIAS$coefficient[5],modelNIAS$coefficient[6],modelNIAS$coefficient[7],modelNIAS$coefficient[8])
pABvec=c(modelNIAS$coefficient[1],modelNIAS$coefficient[2],modelNIAS$coefficient[3],modelNIAS$coefficient[4])



ANIASreport=data.frame(Incentive_Level=incentivec, P_Aresp_in_Astate=pAAvec,P_Aresp_in_Bstate=pABvec,Pval=Pvalvec,correct=fitun)

#Table 6***********************************************
stargazer(ANIASreport,summary=FALSE,rownames=FALSE)
#***********************************************

